home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / prog / yamp2.zip / TEST.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-16  |  6.7 KB  |  219 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/23/92
  7.  
  8. Copyright(c) Mark Von Tress 1992
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission.
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19. // make sure to define the global define variable  IN_RAM if you want
  20. // to use the medium model. This can be done in the compiler code
  21. // generation window of compiler options. It may also be defined for
  22. // the command line compiler
  23.  
  24. #include "virt.h"
  25.  
  26. /* required global declaration for the matrix stack object */
  27.  
  28. //unsigned int _stklen = STACKLENGTH;
  29.  
  30. MStack *Dispatch = new MStack;
  31.  
  32.  
  33. VMatrix& function1( VMatrix &a, VMatrix &b )
  34. {
  35.    // This function tests the freind functions and returns a value
  36.  
  37.    a.Garbage("function1");        // check a and b
  38.    b.Garbage("function1");
  39.  
  40.    Dispatch->Inclevel();          // increment push-pop level
  41.    VMatrix c("c",1,1);            // create a local matrix
  42.  
  43.    c = a + b + a + b;             // check a repeated matrix addition
  44.    c.DisplayMat();                // print c
  45.    c = 431.2 + Tran(a)*b + 2.134; // check commutivity of scalar addition
  46.    c.DisplayMat();                // print c
  47.    c = -c;                        // check uniary minus
  48.    c.DisplayMat();                // print c
  49.    c = a - b;                     // matrix subraction
  50.    c.DisplayMat();                // print c
  51.    c = 5 - a - 5;                 // check commutivity of scalar subtraction
  52.    c.DisplayMat();                // print c
  53.    c = 5*a*5;                     // check commutivity of scalar multiplication
  54.    c.DisplayMat();                // print c
  55.    c = a%a;                       // check elementwise multiplication
  56.    c.DisplayMat();                // print c
  57.    c = a/1234;                    // check scalar division
  58.    c.DisplayMat();                // print c
  59.    c = a/b;                       // check elementwise division
  60.    c.DisplayMat();                // print c
  61.  
  62.    Dispatch->PrintStack();        // show stack before push
  63.    Dispatch->Push(c);             // push c onto stack
  64.    Dispatch->PrintStack();        // examine stack after push
  65.    return Dispatch->DecReturn();  // decrement push-pop level, and return
  66.                   // stack top
  67. }
  68.  
  69. VMatrix &function0 ( void )
  70. {
  71.     // test some of the output functions and raw matrix functions
  72.  
  73.     Dispatch->Inclevel();
  74.     VMatrix d = VMatrix("d",1,1);
  75.     VMatrix a,c,H;
  76.  
  77.     a = Reada( "catchv.dat" );      // read an ascii matrix
  78.     a.DisplayMat();                 // display a
  79.     Writea( "junk.dat", a);         // write ascii matrix
  80.     a.Writeb( "junk.bin", a);       // write binary matrix
  81.     a = Readb( "junk.bin" );        // read binary matrix
  82.     a.InfoMat();                    // display matrix info
  83.     a.DisplayMat();                 // display a
  84.     d = Submat( a, a.r, 4, 1, 2 );  // take a submatrix of a
  85.     d.DisplayMat();                 // display it
  86.     c = Ch( d,d);                   // horizontal concatenation
  87.     c.DisplayMat();
  88.     c = Cv( d,d);                   // vertical contatenation
  89.     c.DisplayMat();
  90.     c = Kron(Ident( 3 ),d);         // Kroniker's product
  91.     Setdec( 1 );                    // set number of decimals to print
  92.     Setwid( 5 );                    // set print width
  93.     c.DisplayMat();
  94.     H = Helm( 4 );                  // make a helmert matrix
  95.     H.DisplayMat();                 // show that it is orthonormal
  96.     d = Tran( H )*H;
  97.     d.DisplayMat();
  98.     d = H*Tran(H);
  99.     d.DisplayMat();
  100.  
  101.     a = Ident( 4 ) + Fill(4,4,0.5); // redefine a
  102.     a.DisplayMat();                 // print a
  103.     c = Tran(Helm(4))*a*H;          // diagonalize a
  104.     c.DisplayMat();                 // display c
  105.     c = Inv( a );                   // invert a
  106.     c.DisplayMat();                 // display c
  107.     VMatrix b;                      // create b
  108.     b = c*a;                        // redefine b
  109.     b.DisplayMat();                 // display b
  110.     c = function1(a,a);             // call a function
  111.     c.DisplayMat();                 // display c
  112.     Dispatch->Push(c);
  113.     return Dispatch->DecReturn();
  114. }
  115.  
  116. VMatrix ®ression( void )  // do a multiple linear regression
  117.   {
  118.       Dispatch->Inclevel();
  119.       VMatrix a,xy,reg;
  120.       a = Reada( "catchv.dat");                 // read data
  121.       VMatrix b = VMatrix("b", a.r, a.c);       //simplify indexes
  122.       int N=a.r;                                // N
  123.       int p=3;                                  // params
  124.       xy = Ch( Fill(N,1,1), Submat(a,N,4,1,2) );// make x y
  125.       reg = Sweep( 1,p, Tran(xy)*xy);           // solve regression
  126.                         // using sweep
  127.       reg.M(p+1,p+1) = reg.m(p+1,p+1)/((double )( N-p));
  128.                         // divide to get mse
  129.  
  130.       VMatrix x= Submat(xy,N,3,1,1);
  131.       VMatrix y= Submat(xy,N,4,1,4);
  132.       VMatrix betahat = Inv(Tran(x)*x)*Tran(x)*y; //solve regression using
  133.                           //normal equations
  134.       betahat.DisplayMat();
  135.       reg.DisplayMat();
  136.       Dispatch->Push( reg );                    // put return mat on stack
  137.       return Dispatch->DecReturn();             // dec level & return
  138. }
  139.  
  140. void altermatrix( VMatrix &t ) // alter a matrix element
  141.   {
  142.      t.M(1,1) = 2525.0;
  143.      VMatrix temp = MSort( t, 1 );
  144.      temp.DisplayMat();
  145.   }
  146.  
  147. void testhuge( void )
  148. {
  149.   VMatrix VeryBig("vb",300,300);
  150.   VeryBig.InfoMat();
  151.  
  152. }
  153.  
  154. void version1p1 ( void )
  155. {
  156.    //
  157.    // testreg.cpp for test of regression functions.
  158.    //
  159.    VMatrix A=Fill(6,6,1.0);
  160.    VMatrix B=Fill(6,1,2.0);
  161.  
  162.    (Vec(A)).DisplayMat();
  163.    (Vecdiag(A)).DisplayMat();
  164.    (Diag(B)).DisplayMat();
  165.    (Shape(A,3)).DisplayMat();
  166.    (Sum(A,ROWS)).DisplayMat();
  167.    (Sumsq(A,COLUMNS)).DisplayMat();
  168.    (Cusum(A)).DisplayMat();
  169.    (Mmin(B)).DisplayMat();
  170.    (Mmax(B,ROWS)).DisplayMat();
  171.  
  172.    VMatrix C = Ch(A,Diag(B));
  173.    Setwid(5);
  174.    Setdec(2);
  175.    C.DisplayMat();
  176.    Crow( C, 1, 0.6 );
  177.    C.DisplayMat();
  178.    Srow( C, 1, 6 );
  179.    C.DisplayMat();
  180.    Lrow( C, 2,3, .5);
  181.    C.DisplayMat();
  182.    Ccol( C, 1, 0.6 );
  183.    C.DisplayMat();
  184.    Scol( C, 1,6 );
  185.    C.DisplayMat();
  186.    Lcol( C, 2,3, .5);
  187.    C.DisplayMat();
  188. }
  189.  
  190. main()
  191. {
  192.   VMatrix testit;
  193.  
  194.   printf("%d\n", sizeof(testit) );
  195.  
  196.   testit = regression();
  197.   testit.DisplayMat();
  198.  
  199.   VMatrix b = testit;
  200.   VMatrix a = 2*testit;
  201.  
  202.   VMatrix c = a + b + a + b;
  203.   c.DisplayMat();
  204.  
  205.   VMatrix t = function0() + function0();
  206.   t.DisplayMat();
  207.  
  208.   altermatrix( testit );
  209.   testit.DisplayMat();
  210.  
  211.   version1p1();
  212.  
  213. #ifndef IN_RAM
  214.        testhuge();
  215. #endif
  216.  
  217.   vclose();
  218. }
  219.